{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Discussion 4\n",
    "\n",
    "In this discussion section we study the learning problem for one-dimensional function $f$ on the segment $[0,1]$ (**which has been studied in Discussion 1**) by using **kernel ridge regression**. It goes in parallel with problem 4 from HW4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.kernel_ridge import KernelRidge\n",
    "from sklearn.metrics.pairwise import pairwise_kernels\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib notebook\n",
    "%matplotlib inline\n",
    "from ipywidgets import interactive\n",
    "import ipywidgets as widgets\n",
    "from ipywidgets import fixed\n",
    "from numpy import linalg as LA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this problem you may encounter warnings that cause jupyter notebook to jump to them abruptly, which is mildly irritating. If you want to suppress those warnings, uncomment the code in the following cell and run it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import warnings\n",
    "# warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#kernel\n",
    "def generate_kernel_widget():\n",
    "    return widgets.Dropdown(\n",
    "        options=['rbf', 'poly'],\n",
    "        description='kernel: ',\n",
    "        disabled=False)\n",
    "\n",
    "#parameter of the kernel\n",
    "def generate_gamma_widget():\n",
    "    return widgets.FloatLogSlider(\n",
    "        value=256,\n",
    "        base=np.sqrt(2),\n",
    "        min=-20,\n",
    "        max=40, \n",
    "        step=1,\n",
    "        description='$\\gamma$',\n",
    "        continuous_update= False)\n",
    "\n",
    "#parameter of poly degree\n",
    "def generate_degree_widget():\n",
    "    return widgets.IntSlider(\n",
    "        value=5,\n",
    "        min=1,\n",
    "        max=100,\n",
    "        step=1,\n",
    "        description='$d$',\n",
    "        disabled=False,\n",
    "        continuous_update=False,\n",
    "        orientation='horizontal',\n",
    "        readout=True,\n",
    "        readout_format='d')\n",
    "\n",
    "#number of samples\n",
    "def generate_n_samples_widget():\n",
    "    return widgets.IntSlider(\n",
    "        value=100,\n",
    "        min=1,\n",
    "        max=100,\n",
    "        step=1,\n",
    "        description='n: ',\n",
    "        disabled=False,\n",
    "        continuous_update=False,\n",
    "        orientation='horizontal',\n",
    "        readout=True,\n",
    "        readout_format='d')\n",
    "    \n",
    "# Selected index\n",
    "def generate_pt_index_widget():\n",
    "    return widgets.IntSlider(\n",
    "        value=50,\n",
    "        max=100-1, \n",
    "        step=1, \n",
    "        description='Index i: ', \n",
    "        disabled=False, \n",
    "        continuous_update=False, \n",
    "        orientation='horizontal', \n",
    "        readout=True, \n",
    "        readout_format='d')\n",
    "\n",
    "# Scale of y\n",
    "def generate_scale_widget():\n",
    "    return widgets.FloatSlider(\n",
    "        value=1.,\n",
    "        min=1.,\n",
    "        max=10,\n",
    "        step=0.01, \n",
    "        description='$y$ scale', \n",
    "        continuous_update= False)\n",
    "\n",
    "#level of regularization\n",
    "def generate_reg_widget():\n",
    "    return widgets.FloatLogSlider(\n",
    "        value=2.11e-8,\n",
    "        base=np.sqrt(2),\n",
    "        min=-60,\n",
    "        max=40, \n",
    "        step=1,\n",
    "        description='Reg ($\\lambda$):',\n",
    "        continuous_update= False)\n",
    "\n",
    "# Noise level\n",
    "def generate_sigma_widget():\n",
    "    return widgets.FloatSlider(\n",
    "        value=0.3,\n",
    "        min=0., \n",
    "        max=0.5, \n",
    "        step=0.01,\n",
    "        description='Noise ($\\sigma$)',\n",
    "        continuous_update= False)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part (a) ##\n",
    "\n",
    "In this problem we study application of kernel ridge regression to the problem of learning a 1d function from samples. We are going to consider two different kernels: rbf and polynomial. The corresponding kernel functions are given by the following:\n",
    "$$\n",
    "K_{\\text{rbf}}(x_i, x_j) = \\exp(-\\gamma(x_i - x_j)^2), \\quad K_{\\text{poly}}(x_i, x_j) = (1+\\gamma x^{\\top}_i x_j)^{d},\n",
    "$$\n",
    "\n",
    "where $\\gamma$ is a parameter of the kernel, and $d$ is a parameter for the degree of poly kernel.\n",
    "\n",
    "**Use the code in the following cell to plot $K(1, x)$ as a function of $x$ for $K$ being 'rbf' or 'poly' kernel. Report how changing parameter $\\gamma$ influences the shape of the plot.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_kernel(kernel, gamma, degree):\n",
    "    if kernel in ['poly']:\n",
    "        kernel_params={'degree': degree,\n",
    "                       'gamma': gamma,\n",
    "                       'coef0': 1.0\n",
    "                      }\n",
    "    else:\n",
    "        kernel_params={'gamma': gamma}\n",
    "    X = np.expand_dims(np.linspace(-5., 5, 1000), axis=1) \n",
    "    k = pairwise_kernels(X, [[1]], kernel, **kernel_params)\n",
    "    plt.xlim(-5, 5)\n",
    "    plt.plot(X[:,0], k)\n",
    "    \n",
    "kernel_choice = generate_kernel_widget()\n",
    "\n",
    "gamma_choice= generate_gamma_widget\n",
    "\n",
    "interactive_plot = interactive(plot_kernel,\n",
    "                               kernel=generate_kernel_widget(),\n",
    "                               gamma=generate_gamma_widget(),\n",
    "                               degree=generate_degree_widget(),     \n",
    "                              )\n",
    "interactive_plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part (b) ##\n",
    "\n",
    "The central object in kernel regression is kernel matrix.  In this part we are going to study its spectral properties. The code in the following cell samples $n$ points from the uniform distribution on $[0,1]$ and plots singular values of the kernel matrix ${\\bf K} =\\{K({\\bf x}_i, {\\bf x}_j)\\}_{i,j=1}^n$. \n",
    "\n",
    "**Before conducting the experiement, report how you think singular values of ${\\bf K}$ will depend on $\\gamma$. Provide some intuition using the results of the previous part. Then run the following cell. Report how singular values of ${\\bf K}$ actually depend on $n$ and $\\gamma$.  Was your initial guess correct?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_eig_values(n_samples, kernel, gamma, degree, seed):\n",
    "    if kernel in ['poly']:\n",
    "        kernel_params={'degree': degree,\n",
    "                   'gamma': gamma,\n",
    "                   'coef0': 1.0\n",
    "                  }\n",
    "    else:\n",
    "        kernel_params={'gamma': gamma}\n",
    "    np.random.seed(seed)\n",
    "    X = np.random.rand(n_samples,1)\n",
    "    K = pairwise_kernels(X, X, kernel, **kernel_params)\n",
    "\n",
    "    eig_vals,_ = np.linalg.eig(K) \n",
    "    eig_vals = np.sort(np.abs(eig_vals))[::-1]\n",
    "\n",
    "    plt.yscale('log')\n",
    "    plt.ylim(1e-20, 1e4)\n",
    "    plt.xlim([-1, n_samples+1])\n",
    "    plt.plot(eig_vals, 'o-', markersize=3.0)\n",
    "    plt.title('eigenvalues of K')\n",
    "             \n",
    "         \n",
    "interactive_plot = interactive(plot_eig_values,\n",
    "                               n_samples=generate_n_samples_widget(),\n",
    "                               kernel=generate_kernel_widget(),\n",
    "                               gamma=generate_gamma_widget(),\n",
    "                               degree=generate_degree_widget(),\n",
    "                               seed=fixed(1))\n",
    "interactive_plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part (c) ##\n",
    "\n",
    "Now we implement the learning procedure. In part a of the problem \"Kernel Ridge Regression: Theory\" we derive that the prediction of ridge regression is given by \n",
    "$$\\hat{f}({\\bf x}) = \\sum_{i=1}^n {\\bf \\hat\\alpha}[i]K({\\bf x} , {\\bf x}_i).$$\n",
    "\n",
    "**Fill in the gap in the code below to compute the weight vector ${\\bf \\hat\\alpha}$**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_predictor(X, y, mode, lambda_reg, kernel, degree=None, gamma=1.0, coef0=1.0):\n",
    "    if mode == 'sklearn':\n",
    "        clf = KernelRidge(alpha=lambda_reg, kernel=kernel, degree=degree, gamma=gamma, coef0=coef0)\n",
    "        clf.fit(X, y)\n",
    "        return clf.predict, clf.dual_coef_\n",
    "    else:\n",
    "        if kernel in ['poly']:\n",
    "            kernel_params={'degree': degree,\n",
    "                           'gamma': gamma,\n",
    "                           'coef0': coef0\n",
    "            }\n",
    "        else:\n",
    "            kernel_params={'gamma': gamma}\n",
    "        \n",
    "        K = pairwise_kernels(X, X, kernel, **kernel_params) + lambda_reg * np.eye(len(y))\n",
    "        ## TODO: solve the dual coefficients 'alpha'\n",
    "        ### start c ###\n",
    "        \n",
    "        ### end c ###\n",
    "        return lambda x: pairwise_kernels(x, X, kernel, **kernel_params) @ alpha, alpha"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part (d) ##\n",
    "\n",
    "We know that our prediction $\\hat{f}(x)$ is linear in ${\\bf y}$ - our measured values. In this part we will see how each individual measurement affects the prediction. \n",
    "\n",
    "The following code allows you to choose index $i$ and then plots the prediction of ridge regression for the case when $y_i = y_{\\text{scale}}$, and $y_{j} = 0$ for $j \\neq i$. In other words, we observe the response that impulse signal in the $i$-th point adds to the total prediction $\\hat{f}({\\bf x}).$ The second plot shows coefficients of the vector $\\hat{\\bf \\alpha}$.\n",
    "\n",
    "1. **Set lambda to the minimal value. For which kernel does the impulse response behave more robustly when you change n and gamma? How does this align with the results of part b?**\n",
    "2. **Report what happens to the vector ${\\bf \\hat\\alpha}$ as you increase gamma. Explain why this happens using part a and the formula for ${\\bf \\hat\\alpha}$**\n",
    "3. **Set lambda to be of order $10^{-7}$. For which value of gamma is the impulse response closest to the step function? How do you think the performance of ridge regression in this regime would compare to 1 nearest neighbour regression?**\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_kernel_impulse(kernel, n_samples, gamma, degree, lambda_reg, pt_index, pt_scale):\n",
    "    n_features  = 1\n",
    "    rng = np.random.RandomState(1)\n",
    "\n",
    "    X = np.sort(rng.rand(n_samples, n_features), axis=0)\n",
    "    \n",
    "    # Impulse (set the [pt_index]-th in y as 1)\n",
    "    y = np.eye(n_samples)[pt_index, :] * pt_scale\n",
    "    \n",
    "    def stair_function(X, value_x):\n",
    "        band_width = 0.02\n",
    "        num = X.shape[0]\n",
    "        Y = np.zeros(num)\n",
    "        idx = np.where(np.abs(X - value_x) < band_width)[0]\n",
    "        Y[idx] = 1.0\n",
    "        return Y\n",
    "\n",
    "    # Kernel Ridge Regression Solver\n",
    "    predictor, alpha_ = compute_predictor(X, y, mode='self-implement', lambda_reg=lambda_reg, kernel=kernel, gamma=gamma, degree=degree)\n",
    "    \n",
    "    clip_bound = 2.5\n",
    "    # Sample test data points\n",
    "    X_test = np.concatenate([X.copy(), np.expand_dims(np.linspace(0., 1., 1000), axis=1)])\n",
    "    X_test = np.sort(X_test, axis=0)\n",
    "\n",
    "    # Visualization\n",
    "    plt.figure(figsize=(10,5))\n",
    "    plt.xlim(0, 1)\n",
    "    plt.ylim(-clip_bound, clip_bound)\n",
    "    plt.plot(X_test, stair_function(X_test, X[pt_index, :]), '--', c='coral', label='stair')\n",
    "    plt.scatter(X[:,0], y, c='blue', s=5, label='Data points')\n",
    "    plt.scatter(X[pt_index,0], y[pt_index], c='darkorange', label='One data point')\n",
    "    plt.plot(X_test, predictor(X_test), '-', color='forestgreen',  label='xi_{i}(x_test)')\n",
    "    plt.title('Impulse function')\n",
    "    plt.legend()\n",
    "    plt.show()\n",
    "    \n",
    "    plt.figure(figsize=(10,6))\n",
    "    plt.xlim(0, 1)\n",
    "    markerlines, stemlines,  baseline = plt.stem(X[:,0], alpha_, 'b', 'o',  use_line_collection=True)\n",
    "    plt.setp(stemlines, 'color', plt.getp(markerlines,'color'))\n",
    "    plt.title('Weights in points')\n",
    "    plt.show()\n",
    "    \n",
    "\n",
    "interactive_plot_impulse = interactive(plot_kernel_impulse,\n",
    "                                       kernel=generate_kernel_widget(),\n",
    "                                       n_samples=generate_n_samples_widget(), \n",
    "                                       gamma=generate_gamma_widget(),\n",
    "                                       degree=generate_degree_widget(),   \n",
    "                                       lambda_reg=generate_reg_widget(),\n",
    "                                       pt_index=generate_pt_index_widget(),\n",
    "                                       pt_scale=generate_scale_widget())\n",
    "interactive_plot_impulse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part (e) ##\n",
    "\n",
    "In this part we study the approximation abilities of our kernels, i.e. we will generate data without noise and see how well the kernel regression can learn. Run the code in the next cell to see three functions that we will try to learn."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f_true(X):\n",
    "    X = X[:, 0]\n",
    "    y = np.zeros(len(X))\n",
    "    y[X < 0.3] = 8 * (X[X < 0.3] + 0.3)**2 - 2.88\n",
    "    y[X > 0.5] = np.sin(2 * np.pi * (X[X> 0.5]-0.5) * 4)\n",
    "    return y\n",
    "\n",
    "X = np.expand_dims(np.linspace(0, 1, 1000), axis=1)\n",
    "\n",
    "plt.figure(figsize=(8,5))\n",
    "plt.plot(X[:, 0], f_true(X))\n",
    "plt.title('f-dis1')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the following cell we will see how well kernel regression can learn the true function. **Before conducting the expreiements report which kernel in your oppinion will be better for each of the funtions above. Use the previous part to provide some intuition supporting your guess. After that, run the following cell and report the following:**\n",
    "1. **Set lambda to the minimum value and try to approximate both functions using both kernels by choosing gamma. For each of the functions report which kernel achieves better approximation.**\n",
    "2. **Now do the same, but set lambda to be of order $10^{-8}$.**\n",
    "3. **Set gamma to be of order $10^2$ and try to approximate both functions using both kernels by choosing lambda. Answer the same question as before.**\n",
    "4. **Report what happens when either $\\gamma$ or $\\lambda$ is chosen to be very large.**\n",
    "5. **Was your initial guess correct?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_kernel_ridge_regression(kernel, n_samples,  gamma, degree, sigma, lambda_reg):\n",
    "    n_features = 1\n",
    "    rng = np.random.RandomState(1)\n",
    "\n",
    "    X = np.sort(rng.rand(n_samples, n_features), axis=0)\n",
    "    \n",
    "    # Generate y\n",
    "    y = f_true(X) + rng.randn(n_samples) * sigma\n",
    "\n",
    "    # Kernel Ridge Regression Solver\n",
    "    predictor, alpha_ = compute_predictor(X, y, mode='self-implement', lambda_reg=lambda_reg, kernel=kernel, gamma=gamma, degree=degree)\n",
    "    \n",
    "    clip_bound = 2.5\n",
    "    # Sample test data points\n",
    "    X_test = np.concatenate([X.copy(), np.expand_dims(np.linspace(0., 1., 1000), axis=1)])\n",
    "    X_test = np.sort(X_test, axis=0)\n",
    "\n",
    "    # Visualization\n",
    "    plt.figure(figsize=(10,7))\n",
    "    plt.xlim(0, 1)\n",
    "    plt.ylim(-clip_bound, clip_bound)\n",
    "    plt.plot(X_test, predictor(X_test), '-', color='forestgreen',  label='prediction')\n",
    "    plt.scatter(X[:,0], y, c='darkorange', s=40.0, label='training data points')\n",
    "    plt.plot(X_test, f_true(X_test), '--', color='royalblue', linewidth=2.0, label='Ground truth')\n",
    "    plt.legend()\n",
    "\n",
    "\n",
    "interactive_plot = interactive(plot_kernel_ridge_regression,\n",
    "                               n_samples=fixed(100),\n",
    "                               kernel=generate_kernel_widget(),\n",
    "                               gamma=generate_gamma_widget(),\n",
    "                               degree=generate_degree_widget(),\n",
    "                               sigma=fixed(0),\n",
    "                               lambda_reg=generate_reg_widget())\n",
    "interactive_plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part (f) ##\n",
    "\n",
    "In this part (and the following parts) we investigate how our predictions are affected by noise. Here we only consider one of two functions from before (namely $f(x) = f_{\\text{dis-1}}(x)$), but we add noise to our data. $f_{\\text{dis-1}}(x)$ is similar to the the function in **Discussion 1**,\n",
    "\n",
    "$$\n",
    "f_{\\text{dis-1}}(x) = \\begin{cases}\n",
    "8(x + 0.3)^2 - 2.88 & x \\in [0, 0.3]\\\\\n",
    "0& x\\in[0.3, 0.5]\\\\\n",
    "4\\sin(2\\pi (x - 0.5))& x \\in [0.5, 1].\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "The data points $\\{x_i\\}_{i=1}^n$ from $[0, 1]$ and learn $f$ from values $\\{f(x_i) + \\epsilon_i\\}$, where $\\{\\epsilon_i\\}_{i=1}^n$ is some additive noise.\n",
    "\n",
    "1. **Set lambda to minimum value and sigma to zero. For each kernel find the range of gamma for which the prediction looks close to $f(x)$. Then set noise to the maximum value and report for which values of gamma the predictions were more persistent to that noise. How do your observations align with the results of part b?**\n",
    "2. **Set noise to maximum and try to find the best prediction for each kernel (you can choose both gamma and lambda). Report for which kernel you could achieve the better result. For which kernel is the prediction more sensitive to the choice of parameters?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_kernel_ridge_regression(kernel, n_samples,  gamma, degree, sigma, f_type, lambda_reg):\n",
    "    n_features = 1\n",
    "    rng = np.random.RandomState(1)\n",
    "\n",
    "\n",
    "\n",
    "    X = np.sort(rng.rand(n_samples, n_features), axis=0)\n",
    "    \n",
    "    # Generate y\n",
    "    y = f_true(X) + rng.randn(n_samples) * sigma\n",
    "\n",
    "    # Kernel Ridge Regression Solver\n",
    "    predictor, alpha_ = compute_predictor(X, y, mode='self-implement', lambda_reg=lambda_reg, kernel=kernel, gamma=gamma, degree=degree)\n",
    "    \n",
    "    clip_bound = 2.5\n",
    "    # Sample test data points\n",
    "    X_test = np.concatenate([X.copy(), np.expand_dims(np.linspace(0., 1., 1000), axis=1)])\n",
    "    X_test = np.sort(X_test, axis=0)\n",
    "\n",
    "    # Visualization\n",
    "    plt.figure(figsize=(10,7))\n",
    "    plt.xlim(0, 1)\n",
    "    plt.ylim(-clip_bound, clip_bound)\n",
    "    plt.plot(X_test, predictor(X_test), '-', color='forestgreen',  label='prediction')\n",
    "    plt.scatter(X[:,0], y, c='darkorange', s=40.0, label='training data points')\n",
    "    plt.plot(X_test, f_true(X_test), '--', color='royalblue', linewidth=2.0, label='Ground truth')\n",
    "    plt.legend()\n",
    "\n",
    "\n",
    "interactive_plot = interactive(plot_kernel_ridge_regression,\n",
    "                               n_samples=fixed(50),\n",
    "                               kernel=generate_kernel_widget(),\n",
    "                               gamma=generate_gamma_widget(),\n",
    "                               degree=generate_degree_widget(),\n",
    "                               sigma=generate_sigma_widget(),\n",
    "                               f_type = fixed('f-dis1'),\n",
    "                               lambda_reg=generate_reg_widget())\n",
    "interactive_plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part (g) ##\n",
    "\n",
    "Now we observe how our parameters influence the generalization error. In this part, **we fix the kernel parameter ($\\gamma=200$)**.\n",
    "\n",
    "We first study the regularization parameter ($\\lambda$).\n",
    "\n",
    "1. **Report what happens to the best regularization parameter ($\\lambda^{\\star}$) when you increase the noise level $\\sigma$ (from $\\sigma=0.0$ to $\\sigma=0.5$).**\n",
    "2. **Was your initial guess correct?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "def calculate_loss(y_pred, y):\n",
    "    \"\"\"Calculate the loss of the model on data X.\"\"\"\n",
    "    n = y_pred.shape[0]\n",
    "    delta_y = y_pred - y\n",
    "    return np.sum(delta_y**2)/n\n",
    "\n",
    "def plot_test_error_lambda(sigma):\n",
    "    # Generate X_train/y_train\n",
    "    train_n_samples, num_features = 100, 1\n",
    "    X_train = np.random.rand(train_n_samples, num_features)\n",
    "    y_train = f_true(X_train) + np.random.randn(train_n_samples) * sigma\n",
    "\n",
    "    # Generate X_test/y_test\n",
    "    test_n_samples = 10000\n",
    "    X_test = np.random.rand(test_n_samples, num_features)\n",
    "    y_test = f_true(X_test)\n",
    "    # Kernel parameter - gamma\n",
    "    gamma = 200.0\n",
    "\n",
    "    # Calculate test loss/norm of dual solution/condition number of (K+\\lambda*I)\n",
    "    loss_test_reg = []\n",
    "    solution_norm_reg = []\n",
    "    # Setup a range of regularization parameters from 1e-9 to 1e5\n",
    "    lambda_reg_set = np.logspace(-10, 5, 50, endpoint=True) \n",
    "    for lambda_reg in lambda_reg_set:\n",
    "        # solve kernel ridge regression\n",
    "        predictor, alpha_= compute_predictor(X_train, y_train, mode='self-implement', lambda_reg=lambda_reg, kernel='rbf', gamma=gamma)\n",
    "        # get prediction on X_test\n",
    "        y_pred = predictor(X_test)\n",
    "        loss_test_reg.append(calculate_loss(y_pred, y_test))\n",
    "        solution_norm_reg.append(LA.norm(alpha_))\n",
    "    # Find best regularization parameter\n",
    "    idx_min = np.argmin(np.array(loss_test_reg))\n",
    "    # Calculate the test loss of all zero predictor\n",
    "    loss_all_zero = calculate_loss(np.zeros_like(y_test), y_test)\n",
    "\n",
    "    # Visualization\n",
    "    plt.figure(figsize=(7,5))\n",
    "    plt.plot(lambda_reg_set, loss_test_reg, \"d:\", color='royalblue', markersize=5.0, label=\"Test Loss\")\n",
    "    plt.axvline(x=lambda_reg_set[idx_min], color='darkorange', linestyle='-', label='Best Regularization')\n",
    "    plt.axhline(y=loss_all_zero, color='red', linestyle='--', linewidth=1.5, label=\"All-zero-solution\")\n",
    "    plt.legend(loc=\"best\")\n",
    "    plt.xscale('log')\n",
    "    plt.yscale('log')\n",
    "    plt.xlabel('Lambda (Regularization)')\n",
    "    plt.ylabel('Test Loss')\n",
    "    plt.title('Test Loss of Different Regularization (Lambda)')\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "    \n",
    "interactive_plot = interactive(plot_test_error_lambda, sigma=generate_sigma_widget())\n",
    "interactive_plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part (h) ##\n",
    "\n",
    "Next we study the kernel parameter ($\\gamma$). In this part, **we fix the regularization parameter ($\\lambda=10^{-8}$)**.\n",
    "\n",
    "1. **Report what happens to the best kernel parameter ($\\gamma^{\\star}$) when you increase the noise level $\\sigma$ (from $\\sigma=0.0$ to $\\sigma=0.5$).**\n",
    "2. **Is there a difference in terms of the shape of the curve when you change the noise level?** **Try to explain what causes the difference.**\n",
    "3. **Was your initial guess correct?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_test_error_gamma(sigma):\n",
    "    # Define ground truth function\n",
    "\n",
    "    # Generate X_train/y_train\n",
    "    train_n_samples, num_features = 100, 1\n",
    "    X_train = np.random.rand(train_n_samples, num_features)\n",
    "    y_train = f_true(X_train) + np.random.randn(train_n_samples) * sigma\n",
    "\n",
    "    # Generate X_test/y_test\n",
    "    test_n_samples = 10000\n",
    "    X_test = np.random.rand(test_n_samples, num_features)\n",
    "    y_test = f_true(X_test)\n",
    "    \n",
    "    # Calculate loss of zero prediction\n",
    "    loss_all_zero = calculate_loss(np.zeros_like(y_test), y_test)\n",
    "    # Setup a range of gamma parameters from 1e-1 to 1e8\n",
    "    gamma_set = np.logspace(0.1, 9, 50, endpoint=True) * 0.1\n",
    "\n",
    "    # Set lambda as 1e-8\n",
    "    lambda_reg = 1e-8\n",
    "\n",
    "    loss_test_gamma = []\n",
    "    condition_number_gamma = []\n",
    "\n",
    "    for gamma in gamma_set:\n",
    "        # solve kernel ridge regression\n",
    "        predictor, _ = compute_predictor(X_train, y_train, mode='self-implement', lambda_reg=lambda_reg, kernel='rbf', gamma=gamma)\n",
    "        # get prediction on X_test\n",
    "        y_pred = predictor(X_test)\n",
    "        loss_test_gamma.append(calculate_loss(y_pred, y_test))\n",
    "    # Find best gamma parameter\n",
    "    idx_min = np.argmin(np.array(loss_test_gamma))\n",
    "\n",
    "    # Visualization\n",
    "    plt.figure(figsize=(8,6))\n",
    "    plt.plot(gamma_set, loss_test_gamma, \"d:\", color='royalblue', markersize=5.0, label=\"Test Loss\")\n",
    "    plt.axvline(x=gamma_set[idx_min], color='darkorange', linestyle='-', label='Best Regularization')\n",
    "    plt.axhline(y=loss_all_zero, color='red', linestyle='--', linewidth=1.5, label=\"All-zero-solution\")\n",
    "    # plt.axhline(y=0.0, color='green', linestyle=':', linewidth=0.8)\n",
    "    plt.legend(loc=\"best\")\n",
    "    plt.xscale('log')\n",
    "    plt.yscale('log')\n",
    "    plt.xlabel('gamma')\n",
    "    plt.ylabel('Test Loss')\n",
    "    plt.title('Test Loss of Different Kernel Parameter (gamma)')\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "interactive_plot = interactive(plot_test_error_gamma,\n",
    "                               sigma=generate_sigma_widget())\n",
    "interactive_plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part (i) ##\n",
    "\n",
    "\n",
    "#### Compare with the ```sklearn``` implementation\n",
    "\n",
    "In this part, we compare our own implementation with the implementations in **```sklearn```** for kernel ridge regression.\n",
    "\n",
    "When the matrix in solving dual problem becomes close to **singular**, ```sklearn``` will using **least-squares solution** instead of performing $(\\mathbf{K}+\\lambda \\mathbf{I})^{-1}$. You could take a look at the implementations in ```sklearn``` [here is the link](https://github.com/scikit-learn/scikit-learn/blob/647fcb1ac13abd8c2eb2554d526f4ad41fee6778/sklearn/linear_model/_ridge.py#L190)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we take a look at the example shown in Part (f) and study the difference between our implementation with the ```sklearn``` implementation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#level of regularization\n",
    "def generate_reg_widget():\n",
    "    return widgets.FloatLogSlider(\n",
    "        value=8.88e-16,\n",
    "        base=np.sqrt(2),\n",
    "        min=-100,\n",
    "        max=40, \n",
    "        step=1,\n",
    "        description='Reg ($\\lambda$):',\n",
    "        continuous_update= False)\n",
    "\n",
    "#parameter of the kernel\n",
    "def generate_gamma_widget():\n",
    "    return widgets.FloatLogSlider(\n",
    "        value=1.00,\n",
    "        base=np.sqrt(2),\n",
    "        min=-20,\n",
    "        max=40, \n",
    "        step=1,\n",
    "        description='$\\gamma$',\n",
    "        continuous_update= False)\n",
    "\n",
    "def plot_kernel_ridge_regression(kernel, n_samples,  gamma, degree, sigma, f_type, lambda_reg):\n",
    "    n_features = 1\n",
    "    rng = np.random.RandomState(1)\n",
    "\n",
    "\n",
    "\n",
    "    X = np.sort(rng.rand(n_samples, n_features), axis=0)\n",
    "    \n",
    "    # Generate y\n",
    "    y = f_true(X) + rng.randn(n_samples) * sigma\n",
    "\n",
    "    # Kernel Ridge Regression Solver\n",
    "    predictor, alpha_ = compute_predictor(X, y, mode='self-implement', lambda_reg=lambda_reg, kernel=kernel, gamma=gamma, degree=degree)\n",
    "\n",
    "    predictor_sklearn, _ = compute_predictor(X, y, mode='sklearn', lambda_reg=lambda_reg, kernel=kernel, gamma=gamma, degree=degree)\n",
    "\n",
    "    \n",
    "    clip_bound = 2.5\n",
    "    # Sample test data points\n",
    "    X_test = np.concatenate([X.copy(), np.expand_dims(np.linspace(0., 1., 1000), axis=1)])\n",
    "    X_test = np.sort(X_test, axis=0)\n",
    "\n",
    "    # Visualization\n",
    "    plt.figure(figsize=(10,7))\n",
    "    plt.xlim(0, 1)\n",
    "    plt.ylim(-clip_bound, clip_bound)\n",
    "    plt.plot(X_test, predictor(X_test), '-', color='forestgreen',  label='prediction')\n",
    "    plt.plot(X_test, predictor_sklearn(X_test), '--', color='coral',  label='prediction (sklearn)')\n",
    "    plt.scatter(X[:,0], y, c='darkorange', s=40.0, label='training data points')\n",
    "    plt.plot(X_test, f_true(X_test), '--', color='royalblue', linewidth=2.0, label='Ground truth')\n",
    "    plt.legend()\n",
    "\n",
    "\n",
    "interactive_plot = interactive(plot_kernel_ridge_regression,\n",
    "                               n_samples=fixed(50),\n",
    "                               kernel=generate_kernel_widget(),\n",
    "                               gamma=generate_gamma_widget(),\n",
    "                               degree=generate_degree_widget(),\n",
    "                               sigma=generate_sigma_widget(),\n",
    "                               f_type = fixed('f-dis1'),\n",
    "                               lambda_reg=generate_reg_widget())\n",
    "interactive_plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we set the default regularization parameter $\\lambda=10^{-15}$, and compare the test loss of our implemented model and ```sklearn``` model when we change the kernel parameter.\n",
    "\n",
    "1. **Report what happens to the test loss when you use the ```sklearn``` to solve the kernel ridge regression problem.**\n",
    "2. **Set ```lambda_reg = 0``` and report what you observe.**\n",
    "3. **Was your initial guess correct?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_test_error_gamma(sigma):\n",
    "    # Define ground truth function\n",
    "\n",
    "    # Generate X_train/y_train\n",
    "    train_n_samples, num_features = 100, 1\n",
    "    X_train = np.random.rand(train_n_samples, num_features)\n",
    "    y_train = f_true(X_train) + np.random.randn(train_n_samples) * sigma\n",
    "\n",
    "    # Generate X_test/y_test\n",
    "    test_n_samples = 10000\n",
    "    X_test = np.random.rand(test_n_samples, num_features)\n",
    "    y_test = f_true(X_test)\n",
    "    \n",
    "    # Calculate loss of zero prediction\n",
    "    loss_all_zero = calculate_loss(np.zeros_like(y_test), y_test)\n",
    "    # Setup a range of gamma parameters from 1e-1 to 1e8\n",
    "    gamma_set = np.logspace(0.1, 9, 50, endpoint=True) * 0.1\n",
    "\n",
    "    # TODO: setup regularization term\n",
    "    # Set lambda as 1e-15\n",
    "    lambda_reg = 1e-15\n",
    "\n",
    "    loss_test_gamma = []\n",
    "    loss_test_gamma_sklearn = []\n",
    "    condition_number_gamma = []\n",
    "\n",
    "    for gamma in gamma_set:\n",
    "        # solve kernel ridge regression\n",
    "        predictor, _ = compute_predictor(X_train, y_train, mode='self-implement', lambda_reg=lambda_reg, kernel='rbf', gamma=gamma)\n",
    "\n",
    "        predictor_sklearn, _ = compute_predictor(X_train, y_train, mode='sklearn', lambda_reg=lambda_reg, kernel='rbf', gamma=gamma)\n",
    "\n",
    "        # get prediction on X_test\n",
    "        y_pred = predictor(X_test)\n",
    "        loss_test_gamma.append(calculate_loss(y_pred, y_test))\n",
    "        y_pred_sklearn = predictor_sklearn(X_test)\n",
    "        loss_test_gamma_sklearn.append(calculate_loss(y_pred_sklearn, y_test))\n",
    "    # Find best gamma parameter\n",
    "    idx_min = np.argmin(np.array(loss_test_gamma))\n",
    "\n",
    "    # Visualization\n",
    "    plt.figure(figsize=(8,6))\n",
    "    plt.plot(gamma_set, loss_test_gamma, \"d:\", color='royalblue', markersize=3.0, label=\"Test Loss\")\n",
    "    plt.plot(gamma_set, loss_test_gamma_sklearn, \"o-\", color='coral', markersize=3.0, label=\"Test Loss\")\n",
    "    plt.axhline(y=loss_all_zero, color='red', linestyle='--', linewidth=1.5, label=\"All-zero-solution\")\n",
    "    plt.legend(loc=\"best\")\n",
    "    plt.xscale('log')\n",
    "    plt.yscale('log')\n",
    "    plt.xlabel('gamma')\n",
    "    plt.ylabel('Test Loss')\n",
    "    plt.title('Test Loss of Different Kernel Parameter (gamma)')\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "interactive_plot = interactive(plot_test_error_gamma,\n",
    "                               sigma=generate_sigma_widget())\n",
    "interactive_plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part (j) ##\n",
    "\n",
    "### $k$-Fold Cross-Validation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on the above two parts, we find that the test loss can be sensitive to hyperparameters. In this part, **we fix the kernel parameter ($\\gamma=200$)** and try to apply $k$-fold cross validation to select the regularization parameter ($\\lambda$). To start with, we use $k=5$. The function that we learn is $f_{\\text{dis-1}}$.\n",
    "\n",
    "Notice here that there is **no noise** in $y$ for the data pair from test dataset, i.e.,\n",
    "\n",
    "$$y_{\\text{test}} = f_{\\text{dis-1}}(x_{\\text{test}}),$$\n",
    "\n",
    "and we have noise for both the training data points and validation data points.\n",
    "\n",
    "1. **Report the the performance of the picked regularization parameter (by $k$-Fold Cross-Validation) on the test data points.**\n",
    "2. **Study the effect of $k$, change $k=5$ to $k=2, 4, 8, 10, 20$, does the picked hyperparameter sensitive to $k$?**\n",
    "2. **Was your initial guess correct?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.RandomState(1)\n",
    "\n",
    "# Set up noise level\n",
    "sigma = 0.2\n",
    "\n",
    "# Generate X_train/y_train\n",
    "train_n_samples, num_features = 100, 1\n",
    "X_train = rng.rand(train_n_samples, num_features)\n",
    "y_train = f_true(X_train) + rng.randn(train_n_samples) * sigma\n",
    "\n",
    "# Generate X_val/y_val (50 data points)\n",
    "val_n_samples = 50\n",
    "X_val = rng.rand(val_n_samples, num_features)\n",
    "y_val = f_true(X_val) + rng.randn(val_n_samples) * sigma\n",
    "\n",
    "# Generate X_test/y_test\n",
    "test_n_samples = 10000\n",
    "X_test = rng.rand(test_n_samples, num_features)\n",
    "y_test = f_true(X_test)\n",
    "\n",
    "# Set gamma = 200\n",
    "gamma = 200\n",
    "\n",
    "loss_all_zero = calculate_loss(np.zeros_like(y_test), y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import KFold\n",
    "\n",
    "# Setup a range of regularization parameters from 1e-5 to 1e5\n",
    "reg_lambda_set = np.logspace(0.1, 10, 50, endpoint=True) * 0.00001\n",
    "\n",
    "# Merge X_train and X_val\n",
    "X_merge = np.concatenate((X_train, X_val), axis=0)\n",
    "y_merge = np.concatenate((y_train, y_val), axis=0)\n",
    "\n",
    "# Setup K-Fold\n",
    "## TODO: Setup you K here\n",
    "K = 5\n",
    "kf = KFold(n_splits=K)\n",
    "kf.get_n_splits(X_merge)\n",
    "loss_k_fold = []\n",
    "\n",
    "for reg_lambda in reg_lambda_set:\n",
    "    loss_avg = 0.0\n",
    "    for train_index, test_index in kf.split(X_merge):\n",
    "        # The current training points are 'X_merge[train_index]'\n",
    "        # The current validation points are 'X_merge[test_index]'\n",
    "\n",
    "        # 1). Estimate the mean with regularization alpha\n",
    "        # 2). Compute the k-fold loss (L_{k-fold}(\\alpha)) \n",
    "        # 3). Add the the k-fold loss (L_{k-fold}(\\alpha)) to the list 'loss_k_fold'\n",
    "        predictor, _ = compute_predictor(X_merge[train_index], y_merge[train_index], mode='self-implement', lambda_reg=reg_lambda, kernel='rbf', gamma=gamma)\n",
    "        # get prediction on X_merge[test_index]\n",
    "        y_pred = predictor(X_merge[test_index])\n",
    "        loss_avg += calculate_loss(y_pred, y_merge[test_index])/K\n",
    "    loss_k_fold.append(loss_avg)\n",
    "\n",
    "# Select the best hyperparameter via k-fold\n",
    "idx_k_fold_min = np.argmin(np.array(loss_k_fold))\n",
    "\n",
    "# Visualization\n",
    "plt.figure(figsize=(8,6))\n",
    "plt.plot(reg_lambda_set, loss_k_fold, \"d:\", color='royalblue', markersize=5.0, label=\"k-fold-Cross-Validation Loss\")\n",
    "plt.axvline(x=reg_lambda_set[idx_k_fold_min], color='darkorange', linestyle='-', label='Selected Regularization')\n",
    "plt.axhline(y=loss_all_zero, color='red', linestyle='--', linewidth=1.5, label=\"All-zero-solution\")\n",
    "plt.legend(loc=\"best\")\n",
    "plt.xscale('log')\n",
    "plt.yscale('log')\n",
    "plt.xlabel('Lambda (Regularization)')\n",
    "plt.ylabel('Validation Loss')\n",
    "plt.title('k-fold-Cross-Validation Loss of Different Regularization (Lambda)')\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we evaluate the models on the test data points. Note that the there is no noise $\\epsilon$ in the test data points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loss_test = []\n",
    "\n",
    "for reg_lambda in reg_lambda_set:\n",
    "    predictor, _ = compute_predictor(X_merge, y_merge, mode='self-implement', lambda_reg=reg_lambda, kernel='rbf', gamma=gamma)\n",
    "    # get prediction on X_test\n",
    "    y_pred = predictor(X_test)\n",
    "    loss_test.append(calculate_loss(y_pred, y_test))\n",
    "\n",
    "# Select the best hyperparameter for test loss\n",
    "idx_test_min = np.argmin(np.array(loss_test))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we visualize the test loss together with the $k$-fold cross-validation loss in the same plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.RandomState(1)\n",
    "\n",
    "# Visualization\n",
    "plt.figure(figsize=(8,6))\n",
    "plt.plot(reg_lambda_set, loss_test, \"o-\", color='forestgreen', markersize=8.0, label=\"Test Loss\")\n",
    "plt.plot(reg_lambda_set, loss_k_fold, \"d:\", color='royalblue', markersize=5.0, label=\"k-fold-Cross-Validation Loss\")\n",
    "plt.axvline(x=reg_lambda_set[idx_k_fold_min], color='darkorange', linestyle='-', label='Selected Regularization')\n",
    "plt.axhline(y=loss_all_zero, color='red', linestyle='--', linewidth=1.5, label=\"All-zero-solution\")\n",
    "plt.legend(loc=\"best\")\n",
    "plt.xscale('log')\n",
    "plt.yscale('log')\n",
    "plt.xlabel('Lambda (Regularization)')\n",
    "plt.ylabel('Loss')\n",
    "plt.title('Compare Test Loss and k-fold-Cross-Validation Loss')\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Congrats! Hope you get a better understanding via this notebook!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
